import numpy as np
import pandas as pd
import xarray as xr
from matplotlib import pyplot as plt
%matplotlib inline
D:\Programs\anaconda3\lib\site-packages\xarray\backends\cfgrib_.py:27: UserWarning: Failed to load cfgrib - most likely there is a problem accessing the ecCodes library. Try `import cfgrib` to get the full error message warnings.warn(
da = xr.open_dataset('NOAA_NCDC_ERSST_v3b_SST.nc',engine='netcdf4')
da
D:\Programs\anaconda3\lib\site-packages\xarray\backends\plugins.py:61: RuntimeWarning: Engine 'cfgrib' loading failed:
Cannot find the ecCodes library
warnings.warn(f"Engine {name!r} loading failed:\n{ex}", RuntimeWarning)
<xarray.Dataset>
Dimensions: (lat: 89, lon: 180, time: 684)
Coordinates:
* lat (lat) float32 -88.0 -86.0 -84.0 -82.0 -80.0 ... 82.0 84.0 86.0 88.0
* lon (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
* time (time) datetime64[ns] 1960-01-15 1960-02-15 ... 2016-12-15
Data variables:
sst (time, lat, lon) float32 ...
Attributes:
Conventions: IRIDL
source: https://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/...
history: extracted and cleaned by Ryan Abernathey for Research Compu...array([-88., -86., -84., -82., -80., -78., -76., -74., -72., -70., -68., -66.,
-64., -62., -60., -58., -56., -54., -52., -50., -48., -46., -44., -42.,
-40., -38., -36., -34., -32., -30., -28., -26., -24., -22., -20., -18.,
-16., -14., -12., -10., -8., -6., -4., -2., 0., 2., 4., 6.,
8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30.,
32., 34., 36., 38., 40., 42., 44., 46., 48., 50., 52., 54.,
56., 58., 60., 62., 64., 66., 68., 70., 72., 74., 76., 78.,
80., 82., 84., 86., 88.], dtype=float32)array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22.,
24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46.,
48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70.,
72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94.,
96., 98., 100., 102., 104., 106., 108., 110., 112., 114., 116., 118.,
120., 122., 124., 126., 128., 130., 132., 134., 136., 138., 140., 142.,
144., 146., 148., 150., 152., 154., 156., 158., 160., 162., 164., 166.,
168., 170., 172., 174., 176., 178., 180., 182., 184., 186., 188., 190.,
192., 194., 196., 198., 200., 202., 204., 206., 208., 210., 212., 214.,
216., 218., 220., 222., 224., 226., 228., 230., 232., 234., 236., 238.,
240., 242., 244., 246., 248., 250., 252., 254., 256., 258., 260., 262.,
264., 266., 268., 270., 272., 274., 276., 278., 280., 282., 284., 286.,
288., 290., 292., 294., 296., 298., 300., 302., 304., 306., 308., 310.,
312., 314., 316., 318., 320., 322., 324., 326., 328., 330., 332., 334.,
336., 338., 340., 342., 344., 346., 348., 350., 352., 354., 356., 358.],
dtype=float32)array(['1960-01-15T00:00:00.000000000', '1960-02-15T00:00:00.000000000',
'1960-03-15T00:00:00.000000000', ..., '2016-10-15T00:00:00.000000000',
'2016-11-15T00:00:00.000000000', '2016-12-15T00:00:00.000000000'],
dtype='datetime64[ns]')[10957680 values with dtype=float32]
weights = np.cos(np.deg2rad(da.lat))
sst_weighted = da.sst.sel(lat=slice(-5,5),lon=slice(190,240)).weighted(weights)
region_mean = sst_weighted.mean(dim=('lon','lat'))
monthly_sst = region_mean.groupby('time.month').mean()
x = region_mean.size
anoma = np.zeros(x-2)
for i in range(0,x-2):
k = i % 12
anoma[i] = region_mean[i] - monthly_sst[k]
## The criteria, that is often used to classify El Niño episodes,
# is that five consecutive 3-month running mean SST anomalies exceed the threshold.
def moving_average(x,w):
return np.convolve(x, np.ones(w),'valid') / w
mean_anomalies = moving_average(anoma,3)
mean_anomalies
# create a dataframe
# 'JFM','FMA','MAM','AMJ','MJJ','JJA','JAS','OND','NDJ','DJF'
alist = ['JFM','FMA','MAM','AMJ','MJJ','JJA','JAS','OND','NDJ','DJF']
templist=[]
for i in range(0,68):
for j in range(0,10):
templist.append(alist[j])
templist.append('JFM')
templist.append('FMA')
df = pd.DataFrame(mean_anomalies,columns=['mean_3_anomalies'])
df['mo'] = pd.DataFrame(templist)
df['date'] = pd.DataFrame(da.time.values)
sub = np.zeros(680)
df
| mean_3_anomalies | mo | date | |
|---|---|---|---|
| 0 | -0.352137 | JFM | 1960-01-15 |
| 1 | -0.307922 | FMA | 1960-02-15 |
| 2 | -0.210943 | MAM | 1960-03-15 |
| 3 | -0.240803 | AMJ | 1960-04-15 |
| 4 | -0.225803 | MJJ | 1960-05-15 |
| ... | ... | ... | ... |
| 675 | 0.500900 | JJA | 2016-04-15 |
| 676 | -0.072014 | JAS | 2016-05-15 |
| 677 | -0.442696 | OND | 2016-06-15 |
| 678 | -0.618626 | NDJ | 2016-07-15 |
| 679 | -0.728378 | DJF | 2016-08-15 |
680 rows × 3 columns
# Visualize the computed Niño 3.4.
df.plot(x="date",y="mean_3_anomalies")
plt.axhline(y=0.5,ls="-",c="red")#添加水平直线
plt.axhline(y=0,ls=":",c="black")#添加水平直线
plt.axhline(y=-0.5,ls="-",c="blue")#添加水平直线
plt.xticks(rotation=45)
plt.ylabel('sst_anomalies[℃]')
plt.xlabel('year')
plt.title('Niño 3.4(5°N-5°S,120°W-170°W)')
plt.show()
d = xr.open_dataset('CERES_EBAF-TOA_200003-201701.nc',engine='netcdf4')
fig = plt.figure(figsize=(10,6), dpi=120)
grid = plt.GridSpec(5, 9) # 3 rows 3 cols
plt.subplot(grid[0:2, 0:3])
ax1 = d.toa_sw_all_mon.mean(dim='time').plot()
plt.title('TOA Shortwave Flux - All-Sky',fontsize=10)
plt.subplot(grid[0:2, 3:6])
ax2 = d.toa_lw_all_mon.mean(dim='time').plot()
plt.title('TOA Longwave Flux - All-Sky',fontsize=10)
plt.subplot(grid[0:2, 6:9])
ax3 = d.solar_mon.mean(dim='time').plot()
plt.title('Incoming Solar Flux',fontsize=10)
plt.subplot(grid[2:5, 0:5])
ax4 = d.toa_net_all_mon.mean(dim='time').plot()
plt.title('TOA Net Flux - All-Sky')
plt.subplot(grid[2:5, 4:10])
ax5 = (d.solar_mon.mean(dim='time') - d.toa_sw_all_mon.mean(dim='time') - d.toa_lw_all_mon.mean(dim='time') ).plot()
plt.title(' Add up three variables')
plt.tight_layout()
[Hint: Consider calculating the area of each grid]
temp=[]
x = d.solar_mon.mean(dim='time') - d.toa_sw_all_mon.mean(dim='time') - d.toa_lw_all_mon.mean(dim='time')
y = d.toa_net_all_mon.mean(dim='time')
a = plt.scatter(x,y)
plt.xlabel('calculate_toa',fontsize=15)
plt.ylabel('toa_net',fontsize=15)
(x-y).sum()
<xarray.DataArray ()> array(1.6990367, dtype=float32)
array(1.6990367, dtype=float32)
weights = np.cos(np.deg2rad(d.lat))
net_weighted = d.toa_net_all_mon.weighted(weights)
region_mean = net_weighted.sum(dim=['lon','time']).plot()
plt.xlabel('lat/°',fontsize=15)
plt.ylabel('toa_net_all_mon/ W',fontsize=15)
plt.title('Total amount of net radiation')
plt.show()
fig = plt.figure(figsize=(10,6), dpi=120)
grid = plt.GridSpec(4, 4) # 4 rows 4 cols
plt.subplot(grid[0:2, 0:2])
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon<=25)).mean(dim='time').plot()
plt.title('sw-low cloud area regions')
plt.subplot(grid[0:2, 2:4])
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon>=75)).mean(dim='time').plot()
plt.title('sw-high cloud area regions')
plt.subplot(grid[2:4, 0:2])
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon<=25)).mean(dim='time').plot()
plt.title('lw-low cloud area regions')
plt.subplot(grid[2:4, 2:4])
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon>=75)).mean(dim='time').plot()
plt.title('lw-high cloud area regions')
plt.tight_layout()
fig = plt.figure(figsize=(18,6), dpi=600)
grid = plt.GridSpec(1, 2) # 1 rows 2 cols
plt.subplot(grid[0, 0])
d.toa_sw_all_mon.weighted(weights).mean(dim=['lon','lat']).plot(label='global')
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon>=75)).weighted(weights).mean(dim=['lon','lat']).plot(label='high cloud region')
d.toa_sw_all_mon.where((d.cldarea_total_daynight_mon<=25)).weighted(weights).mean(dim=['lon','lat']).plot(label='low cloud region')
plt.legend(loc='best')
plt.xlabel('time',fontsize=18)
plt.ylabel('toa_sw_all_mon',fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.subplot(grid[0, 1])
d.toa_lw_all_mon.weighted(weights).mean(dim=['lon','lat']).plot(label='global')
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon>=75)).weighted(weights).mean(dim=['lon','lat']).plot(label='high cloud region')
d.toa_lw_all_mon.where((d.cldarea_total_daynight_mon<=25)).weighted(weights).mean(dim=['lon','lat']).plot(label='low cloud region')
plt.legend(loc='best')
plt.xlabel('time',fontsize=18)
plt.ylabel('toa_lw_all_mon',fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.tight_layout()
T = xr.open_dataset('2016-2020-T.nc',engine="netcdf4")
T.t2m.mean(dim=('latitude', 'longitude'))
<xarray.DataArray 't2m' (time: 43848)>
array([286.67023, 288.65994, 289.85648, ..., 278.07663, 277.76575,
277.50546], dtype=float32)
Coordinates:
* time (time) datetime64[ns] 2016-01-01 ... 2020-12-31T23:00:00array([286.67023, 288.65994, 289.85648, ..., 278.07663, 277.76575,
277.50546], dtype=float32)array(['2016-01-01T00:00:00.000000000', '2016-01-01T01:00:00.000000000',
'2016-01-01T02:00:00.000000000', ..., '2020-12-31T21:00:00.000000000',
'2020-12-31T22:00:00.000000000', '2020-12-31T23:00:00.000000000'],
dtype='datetime64[ns]')weight = np.cos(np.deg2rad(T.latitude))
T_weighted = T.t2m.weighted(weight)
region_mean = T_weighted.mean(dim=('longitude','latitude'))
a = region_mean.values
index = pd.date_range('2016-01-01', periods=1827, freq='D')
k = 0
t_day = np.zeros(1827)
for i in range(0,1827):
t_day[i]=a[k:k+24].mean()
k+=24
tempdf = pd.DataFrame(t_day,columns=['t2m'])
tempdf['date'] = index
tempdf['year'] = pd.DatetimeIndex(tempdf.date).year
tempdf['month'] = pd.DatetimeIndex(tempdf.date).month
t_month = tempdf.groupby(['year','month']).mean()
monthly_T = region_mean.groupby('time.month').mean()
sea = np.zeros(60)
for i in range(0,5):
sea[12*i:12*i+12] = monthly_T
certain = t_month['t2m'] - sea
certain.plot()
plt.ylabel('temperature[℃]')
Text(0, 0.5, 'temperature[℃]')
T.t2m.mean(dim='time').plot()
<matplotlib.collections.QuadMesh at 0x1fcf932bc10>
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from scipy.stats import norm
np.random.seed(0)
# example data
mu = 100 # mean of distribution
sigma = 15 # standard deviation of distribution
x = a
num_bins = 50
fig, ax = plt.subplots()
# the histogram of the data
n, bins, patches = ax.hist(x, num_bins, density=1)
# Tweak spacing to prevent clipping of ylabel
plt.title('Histogram of T: $\mu=100$, $\sigma=15$')
plt.xlabel('temperature')
plt.ylabel('probability dendity')
fig.tight_layout()
plt.show()
data = a
# basic plot
plt.boxplot(data)
{'whiskers': [<matplotlib.lines.Line2D at 0x1fc8d7f9610>,
<matplotlib.lines.Line2D at 0x1fc8d7f99d0>],
'caps': [<matplotlib.lines.Line2D at 0x1fc8d7f9d60>,
<matplotlib.lines.Line2D at 0x1fc8d807130>],
'boxes': [<matplotlib.lines.Line2D at 0x1fc8d7f93d0>],
'medians': [<matplotlib.lines.Line2D at 0x1fc8d8074c0>],
'fliers': [<matplotlib.lines.Line2D at 0x1fc8d807850>],
'means': []}
T_weighted.mean(dim=('longitude','latitude')).plot()
plt.ylabel('t2m[K]')
Text(0, 0.5, 't2m[K]')
monthly_T.plot()
plt.ylabel('t2m[K]')
Text(0, 0.5, 't2m[K]')